home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 26 / Cream of the Crop 26.iso / os2 / octa209s.zip / octave-2.09 / libcruft / linpack / zgedi.f < prev    next >
Text File  |  1996-07-19  |  4KB  |  136 lines
  1.       subroutine zgedi(a,lda,n,ipvt,det,work,job)
  2.       integer lda,n,ipvt(1),job
  3.       complex*16 a(lda,1),det(2),work(1)
  4. c
  5. c     zgedi computes the determinant and inverse of a matrix
  6. c     using the factors computed by zgeco or zgefa.
  7. c
  8. c     on entry
  9. c
  10. c        a       complex*16(lda, n)
  11. c                the output from zgeco or zgefa.
  12. c
  13. c        lda     integer
  14. c                the leading dimension of the array  a .
  15. c
  16. c        n       integer
  17. c                the order of the matrix  a .
  18. c
  19. c        ipvt    integer(n)
  20. c                the pivot vector from zgeco or zgefa.
  21. c
  22. c        work    complex*16(n)
  23. c                work vector.  contents destroyed.
  24. c
  25. c        job     integer
  26. c                = 11   both determinant and inverse.
  27. c                = 01   inverse only.
  28. c                = 10   determinant only.
  29. c
  30. c     on return
  31. c
  32. c        a       inverse of original matrix if requested.
  33. c                otherwise unchanged.
  34. c
  35. c        det     complex*16(2)
  36. c                determinant of original matrix if requested.
  37. c                otherwise not referenced.
  38. c                determinant = det(1) * 10.0**det(2)
  39. c                with  1.0 .le. cabs1(det(1)) .lt. 10.0
  40. c                or  det(1) .eq. 0.0 .
  41. c
  42. c     error condition
  43. c
  44. c        a division by zero will occur if the input factor contains
  45. c        a zero on the diagonal and the inverse is requested.
  46. c        it will not occur if the subroutines are called correctly
  47. c        and if zgeco has set rcond .gt. 0.0 or zgefa has set
  48. c        info .eq. 0 .
  49. c
  50. c     linpack. this version dated 08/14/78 .
  51. c     cleve moler, university of new mexico, argonne national lab.
  52. c
  53. c     subroutines and functions
  54. c
  55. c     blas zaxpy,zscal,zswap
  56. c     fortran dabs,dcmplx,mod
  57. c
  58. c     internal variables
  59. c
  60.       complex*16 t
  61.       double precision ten
  62.       integer i,j,k,kb,kp1,l,nm1
  63. c
  64.       complex*16 zdum
  65.       double precision cabs1
  66.       double precision dreal,dimag
  67.       complex*16 zdumr,zdumi
  68.       dreal(zdumr) = zdumr
  69.       dimag(zdumi) = (0.0d0,-1.0d0)*zdumi
  70.       cabs1(zdum) = dabs(dreal(zdum)) + dabs(dimag(zdum))
  71. c
  72. c     compute determinant
  73. c
  74.       if (job/10 .eq. 0) go to 70
  75.          det(1) = (1.0d0,0.0d0)
  76.          det(2) = (0.0d0,0.0d0)
  77.          ten = 10.0d0
  78.          do 50 i = 1, n
  79.             if (ipvt(i) .ne. i) det(1) = -det(1)
  80.             det(1) = a(i,i)*det(1)
  81. c        ...exit
  82.             if (cabs1(det(1)) .eq. 0.0d0) go to 60
  83.    10       if (cabs1(det(1)) .ge. 1.0d0) go to 20
  84.                det(1) = dcmplx(ten,0.0d0)*det(1)
  85.                det(2) = det(2) - (1.0d0,0.0d0)
  86.             go to 10
  87.    20       continue
  88.    30       if (cabs1(det(1)) .lt. ten) go to 40
  89.                det(1) = det(1)/dcmplx(ten,0.0d0)
  90.                det(2) = det(2) + (1.0d0,0.0d0)
  91.             go to 30
  92.    40       continue
  93.    50    continue
  94.    60    continue
  95.    70 continue
  96. c
  97. c     compute inverse(u)
  98. c
  99.       if (mod(job,10) .eq. 0) go to 150
  100.          do 100 k = 1, n
  101.             a(k,k) = (1.0d0,0.0d0)/a(k,k)
  102.             t = -a(k,k)
  103.             call zscal(k-1,t,a(1,k),1)
  104.             kp1 = k + 1
  105.             if (n .lt. kp1) go to 90
  106.             do 80 j = kp1, n
  107.                t = a(k,j)
  108.                a(k,j) = (0.0d0,0.0d0)
  109.                call zaxpy(k,t,a(1,k),1,a(1,j),1)
  110.    80       continue
  111.    90       continue
  112.   100    continue
  113. c
  114. c        form inverse(u)*inverse(l)
  115. c
  116.          nm1 = n - 1
  117.          if (nm1 .lt. 1) go to 140
  118.          do 130 kb = 1, nm1
  119.             k = n - kb
  120.             kp1 = k + 1
  121.             do 110 i = kp1, n
  122.                work(i) = a(i,k)
  123.                a(i,k) = (0.0d0,0.0d0)
  124.   110       continue
  125.             do 120 j = kp1, n
  126.                t = work(j)
  127.                call zaxpy(n,t,a(1,j),1,a(1,k),1)
  128.   120       continue
  129.             l = ipvt(k)
  130.             if (l .ne. k) call zswap(n,a(1,k),1,a(1,l),1)
  131.   130    continue
  132.   140    continue
  133.   150 continue
  134.       return
  135.       end
  136.